! subroutine dgesv (integer N, integer NRHS, double precision, dimension(
!     lda, * ) A, integer LDA, integer, dimension( * ) IPIV, double
!    precision, dimension( ldb, * ) B, integer LDB, integer INFO)
!    DGESV computes the solution to system of linear equations A * X = B
!  for GE matrices

module interface_d
    interface
        subroutine read_matrix(A, B)
            implicit none
            integer, parameter :: dbl = 8
            real(dbl), allocatable, dimension(:, :), intent(inout) :: A
            real(dbl), allocatable, dimension(:, :), intent(inout) :: B
        end subroutine read_matrix

        subroutine print_matrix(A)
            integer, parameter :: dbl = 8
            real(dbl), dimension(:, :), intent(in) :: A
        end subroutine print_matrix
    end interface
end module interface_d

program dgesv_

    use interface_d
    implicit none
    integer, parameter :: dbl = 8, fnum = 4
    real(dbl), allocatable, dimension(:, :) :: A, B
    integer :: n, nrhs, lda, ldb, info
    integer, allocatable, dimension(:) :: ipiv
    call read_matrix(A, B)

    call print_matrix(A)
    call print_matrix(B)

    n = size(A, 1)
    nrhs = size(B, 2)
    ldb = n; lda = n
    allocate (ipiv(n))
    call dgesv(n, nrhs, A, lda, ipiv, B, ldb, info)
    if (info == 0) then
        write (*, *) 'ok'
    end if
    call print_matrix(B)
    deallocate (A)
    deallocate (B)
    deallocate (ipiv)
end program dgesv_

subroutine read_matrix(A, B)
    implicit none

    integer, parameter :: dbl = 8, fnum = 4
    integer :: dim, b_col
    real(dbl), allocatable, dimension(:, :), intent(inout) :: A
    real(dbl), allocatable, dimension(:, :), intent(inout) :: B
    integer :: k, status

    open (unit=fnum, file='matrix.txt', status='old', action='read', iostat=status)
    openif: if (status /= 0) then
        write (*, *) 'open file error'
        stop
    else
        read (fnum, *) dim, b_col
        write (*, *) dim, b_col
        allocate (A(dim, dim))
        allocate (B(dim, b_col))

        readA: do k = 1, dim
            read (fnum, *) A(k, :)
        end do readA

        readB: do k = 1, b_col
            read (fnum, *) B(:, k)
        end do readB
    end if openif

    close (fnum)
end subroutine read_matrix

! 按行打印
subroutine print_matrix(A)
    integer, parameter :: dbl = 8
    real(dbl), dimension(:, :), intent(in) :: A
    integer :: row, k
    row = size(A, 1)
    do k = 1, row
        write (*, *) A(k, :)
    end do
    write (*, *)
end subroutine print_matrix
